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ABSTRACT 

The dynamical reaction of the particles accelerated at a shock front by the first order 
Fermi process can be determined within kinetic models that account for both the 
hydrodynamics of the shocked fluid and the transport of the accelerated particles. 
These models predict the appearance of multiple solutions, all physically allowed. We 
discuss here the role of injection in selecting the real solution, in the framework of a 
simple phenomenological recipe, which is a variation of what is sometimes referred to 
as thermal leakage. In this context we show that multiple solutions basically disappear 
and when they are present they are limited to rather peculiar values of the parameters. 
We also provide a quantitative calculation of the efficiency of particle acceleration at 
cosmic ray modified shocks and we identify the fraction of energy which is advected 
downstream and that of particles escaping the system from upstream infinity at the 
maximum momentum. The consequences of efficient particle acceleration for shock 
heating are also discussed. 
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1 INTRODUCTION 



Diffusive particle acceleration at non-relativistic shock fronts 
is an extensively studied phenomenon. Detailed discussions 
of the current status of the investigations can be found in 
some excellent reviews (Drury 1983; Blandford & Eichler 
1987; Berezhko & Krimsky 1988; Jones & Ellison 1991; 
Malkov & Drury 2001). While much is by now well un- 
derstood, some issues are still subjects of much debate, for 
the theoretical and phenomenological implications that they 
may have. One of the most important of these is the reac- 
tion of the accelerated particles on the shock: the violation 
of the test particle approximation occurs when the accelera- 
tion process becomes sufficiently efficient that the pressure 
of the accelerated particles is comparable with the incoming 
gas kinetic pressure. Both the spectrum of the particles and 
the structure of the shock are changed by this phenomenon, 
which is therefore intrinsically nonlinear. 

At present there are three viable approaches to deter- 
mine the reaction of the particles upon the shock: one is 
based on the ever-improving numerical simulations (Jones 
& Ellison 1991; Bell 1987; Ellison, Mobius & Paschmann 
1990; Ellison, Baring & Jones 1995, 1996; Kang & Jones 
1997; Kang, Jones & Gieseler 2002; Kang & Jones 2005) 
that allow one to achieve a self-consistent treatment of sev- 
eral effects. 

The second approach is based on the so-called two-fluid 



model, and treats cosmic rays as a relativistic second fluid. 
This class of models was proposed and discussed in (Drury 
& Volk 1980, 1981; Drury, Axford & Summers 1982; Ax- 
ford, Leer & McKenzie 1982; Duffy, Drury & Volk 1994). 
These models allow one to obtain the thermodynamics of 
the modified shocks, but do not provide information about 
the spectrum of accelerated particles. 

The third approach is semi-analytical and may be very 
helpful to understand the physics of the nonlinear effects in 
a way that sometimes is difficult to achieve through sim- 
ulations, due to their intrinsic complexity and limitations 
in including very different spatial scales. Blandford (1980) 
proposed a perturbative approach in which the pressure of 
accelerated particles was treated as a small perturbation. 
By construction this method provides the description of the 
reaction only for weakly modified shocks. 

Alternative approaches were proposed by Eichler 
(1984), Ellison & Eichler (1984), Eichler (1985) and Ellison 
& Eichler (1985), based on the assumption that the diffusion 
of the particles is sufficiently energy dependent that differ- 
ent parts of the fluid are affected by particles with different 
energies. The way the calculations are carried out implies a 
sort of separate solution of the transport equation for sub- 
relativistic and relativistic particles, so that the two spectra 
must be somehow connected at p ~ mc a posteriori. 

In (Berezhko, Yelshin & Ksenofontov 1994; Berezhko, 
Ksenofontov & Yelshin 1995; Berezhko 1996) the effects of 
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the non-linear reaction of accelerated particles on the maxi- 
mum achievable energy were investigated, together with the 
effects of geometry. The maximum energy of the particles 
accelerated in supernova remnants in the presence of large 
acceleration efficiencies was also studied by Ptuskin & Zi- 
rakashvili (2003a,b). 

The need for a practical solution of the accelera- 
tion problem in the non-linear regime was recognized by 
Berezhko & Ellison (1999), where a simple analytical broken- 
power-law approximation of the non-linear spectra was pre- 
sented. 

Recently, some promising analytical solutions of the 
problem of non-linear shock acceleration have appeared in 
the literature (Malkov 1997; Malkov, Diamond & Volk 2000; 
Blasi 2002, 2004). Blasi (2004) considered for the first time 
the important effect of seed pre-existing particles in the ac- 
celeration region (the linear theory of this phenomenon was 
first studied by Bell (1978)). In a recent work by Kang & 
Jones (2005) the seed particles were included in numerical 
simulations of the acceleration process. 

Numerical simulations have been instrumental to iden- 
tify the dramatic effects of the particles reaction: they 
showed that even when the fraction of particles injected 
from the thermal gas is relatively small, the energy chan- 
nelled into these few particles can be an appreciable part of 
the kinetic energy of the unshocked fluid, making the test 
particle approach unsuitable. The most visible effects of the 
reaction of the accelerated particles on the shock appear 
in the spectrum of the accelerated particles, which shows a 
peculiar hardening at the highest energies. The analytical 
approaches reproduce well the basic features arising from 
nonlinear effects in shock acceleration. 

There is an important point which is still lacking in the 
calculations of the non-linear particle acceleration at shock 
waves, namely the possible amplification of the background 
magnetic field, found in the numerical simulations by Lucek 
& Bell (2000, 2000a) and Bell & Lucek (2001) and recently 
described by Bell (2004) . This effect is still ignored in all cal- 
culations of the reaction of cosmic rays on the shock struc- 
ture. We will not include this effect in the present paper. 

Nonlinear effects in shock acceleration of thermal parti- 
cles result in the appearance of multiple solutions in certain 
regions of the parameter space. This phenomenon is very 
general and was found in both the two-fluid models (Drury 
& Volk 1980, 1981) and in the kinetic models (Malkov 1997; 
Malkov et al. 2001; Blasi 2004). Monte Carlo approaches do 
not show multiple solutions. 

This behaviour resembles that of critical systems, with 
a bifurcation occurring when some threshold is reached in 
a given order parameter. In the case of shock acceleration, 
it is not easy to find a way of discriminating among the 
multiple solutions when they appear. In (Mond & Drury 
1998), a two fluid approach has been used to demonstrate 
that when three solutions appear, the one with intermedi- 
ate efficiency for particle acceleration is unstable to corruga- 
tions in the shock structure and emission of acoustic waves. 
Plausibility arguments may be put forward to justify that 
the system made of the shock plus the accelerated particles 
may sit at the critical point (see for instance the paper by 
Malkov, Diamond & Volk (2000)), but we are not aware of 
any real proof that this is what happens. The physical pa- 
rameters that play a role in this approach to criticality are 



the maximum momentum achievable by the particles in the 
acceleration process, the Mach number of the shock, and the 
injection efficiency, namely the fraction of thermal particles 
crossing the shock that are accelerated to nonthermal ener- 
gies. The last of them, the injection efficiency, hides a crucial 
physics problem by itself, and plays an important role in es- 
tablishing the level of shock modification. This efficiency pa- 
rameter in reality is determined by the microphysics of the 
shock and should not be a free parameter of the problem. 
Unfortunately, our poor knowledge of such microphysics, in 
particular for collisionless shocks, does not allow us to estab- 
lish a clear and universal connection between the injection 
efficiency and the macroscopic shock properties. Put aside 
the possibility to have a fully self-consistent picture of this 
phenomenon, one can try to achieve a phenomenological de- 
scription of it. Kang, Jones & Gieseler (2002) introduced 
a sort of weight function to determine a return probabil- 
ity of particles in the downstream fluid to the upstream 
fluid, as a function of particle momentum. Only sufficiently 
suprathermal particles can jump back to the upstream re- 
gion and therefore take part in the acceleration process. Here 
we adopt an injection recipe which is similar to the thermal 
leakage model of Kang et al. (2002) (see also previous papers 
by Malkov (1998) and by Gieseler et al. (2000)) and imple- 
ment it in the semi-analytical approach of Blasi (2002, 2004). 
We investigate then the phenomenon of multiple solutions 
and show that the injection recipe dramatically reduces the 
appearance of these situations. We also study in some detail 
the efficiency for particle acceleration as a function of the 
Mach number of the shock and the maximum momentum of 
the accelerated particles. 

The paper is structured as follows: in section 2 we briefly 
describe the method proposed by Blasi (2002) for the calcu- 
lation of the spectum and pressure of particles accelerated 
at a modified shock. We describe the appearance of multiple 
solutions in section 3, and the comparison with the method 
of Malkov (1997) in section 4. In section 5 we introduce a 
recipe for the injection of particles from the thermal pool. 
This recipe is then used in section 6 to show how the regions 
of parameter space where multiple solutions appear are dras- 
tically reduced by the self-regulated injection. In section 7 
we discuss the efficiency of particle acceleration at modi- 
fied shocks, and stress the role of escape of particles from 
upstream infinity. The consequences of the cosmic ray mod- 
ification on the shock heating are investigated in section 8. 
We conclude in section 9. 



2 A KINETIC MODEL FOR PARTICLE 

ACCELERATION AT MODIFIED SHOCKS 

In this section we describe the method proposed by Blasi 
(2002, 2004) for the calculation of the spectrum and pres- 
sure of the particles accelerated at a shock surface, when 
the reaction of the particles is taken into account. No seed 
particles are included here. 

The equation that describes the diffusive transport of 
particles in one dimension is 

df(x,p) 



dx 



D m f{x,p) 



dx 



- + 



1 du df(x,p) 



3 dx 



V 



dp 



+ Q(x,p) =0, 



(1) 
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Figure 1. Schematic view of the shock region. 

where we assumed stationarity (df/dt — 0). The x axis is 
oriented from upstream to downstream, as in fig. 1. The 
pressure of the accelerated particles slows down the fluid 
upstream before it crosses the shock surface, therefore the 
gas velocity at upstream infinity, uo, is different from ui, the 
fluid speed just upstream of the shock. The injection term 
is taken in the form Q(x,p) = Qo(p)S(x). 

As a first step, we integrate Eq. 1 around x = 0, from 
x — CP to x = + , denoted in fig. 1 as points "1" and "2" 
respectively, so that the following equation can be written: 









. dx_ 


2 


. dx . 



, 1 df 
i 3 dp 



• ui) + Qo(p) = 0, (2) 



where u\ (U2) is the fluid speed immediately upstream 
(downstream) of the shock and fo is the particle distribution 
function at the shock location. By requiring that the distri- 
bution function downstream is independent of the spatial 
coordinate (homogeneity), we obtain [-D§^] 2 = 0, so that 
the boundary condition at the shock can be rewritten as 



^ dx 



(3) 



We can now perform the integration of Eq. 1 from x = — oo 
to x — CP (point "1"), in order to take into account the 
boundary condition at upstream infinity. Using Eq. 3 we 
obtain 



1 df 
3 P ^ (M2 



Ul) - Ml/o + Qo(p) + 



f 



i p du 

dxf Tx + 



If 



du d f 

dx l~ p JT = °- 
dx op 



We introduce the quantity u p defined as 



1 f° 

U p = Ml - — / 

To J_, 



dx^f(x,p), 



(4) 



(5) 



whose physical meaning is instrumental to understand the 
nonlinear reaction of particles. The function u p is the aver- 
age fluid velocity experienced by particles with momentum 
p while diffusing upstream away from the shock surface. In 
other words, the effect of the average is that, instead of a 
constant speed «i upstream, a particle with momentum p 
experiences a spatially variable speed, due to the pressure of 



the accelerated particles. Since the diffusion coefficient is in 
general p-dependent, particles with different energies feel a 
different compression coefficient, higher at higher energies if, 
as expected, the diffusion coefficient is an increasing func- 
tion of momentum (see (Blasi 2004) for further details on 
the meaning of the quantity u p ). 

With the introduction of u p , Eq. 4 becomes 



1 df 

_ p _ (M2 _ Wp )_ /o 



1 du p 
Up+ 3 P lip- 



+ Qo(p) = o, 



(6) 



where we used the fact that 

Ml 



V 



du 

dx — f = p 

ax 



dfo, . , du p 

^(«i-«p)-/o-^- 



The solution of Eq. 6 for a monochromatic injection at mo- 
mentum Pi n j is 



/o(p) = 



/ 

J Pi 



p dp 3Qo(p) 

p Up — U2 



■ exp 



p' U pl - U 2 



Up, + -p 



, dUp, 
lip 1 " 



3i? s ub f]^gas,l 



R sub - 1 471-pf 



x exp 



■ U 2 



+ 3 P 



1 , dUp, 



dp' 



(7) 



1"gm,l' 



Anp z 



-6{p — Pinj), with n gas ,i the 



Here we used Qo(p) = 

gas density immediately upstream (x — _ ) and 7/ the frac- 
tion of the particles crossing the shock which take part in 
the acceleration process. 

Here we introduced the two quantities R au b = 1*1/1*2 
and Rtot = U0/U2, which are respectively the compression 
factor at the gas subshock and the total compression fac- 
tor between upstream infinity and downstream. For a mod- 
ified shock, Rtot can attain values much larger than R su b 
and more in general, much larger than 4, which is the maxi- 
mum value achievable for an ordinary strong non-relativistic 
shock. The increase of the total compression factor com- 
pared with the prediction for an ordinary shock is respon- 
sible for the peculiar flattening of the spectra of acceler- 
ated particles that represents a feature of nonlinear effects 
in shock acceleration. In terms of R su b and R to t the density 
immediately upstream is n gas ,i = (po/m p )R t ot/ 'Rsub- 

In Eq. 7 we can introduce a dimensionless quantity 
U(p) — Up/uo so that 



x exp ■ 



ft j (p) 



p' 



3R S 



3R 



RtotU(p) 
U(p') 



RtotU{ P >) - 1 



(8) 



The nonlinearity of the problem reflects in the fact that 
U (p) is in turn a function of fo as it is clear from the defini- 
tion of u p . In order to solve the problem we need to write the 
equations for the thermodynamics of the system including 
the gas, the cosmic rays accelerated from the thermal pool 
and the shock itself. 

The velocity, density and thermodynamic properties of 
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the fluid can be determined by the mass and momentum 
conservation equations, with the inclusion of the pressure of 
the accelerated particles. We write these equations between 
a point far upstream (x = — oo), where the fluid velocity is 
Mo and the density is po = mn gas fi, and the point where the 
fluid velocity is u v (density p p ). The index p denotes quan- 
tities measured at the point where the fluid velocity is u p , 
namely at the point x p that can be reached only by parti- 
cles with momentum > p (this is clearly an approximation, 
but as shown in section 4 it provides a good agreement with 
other calculations where this approximation is not used). 
The mass conservation implies: 

pouo = p P u p . (9) 

Conservation of momentum reads: 

POUq + Pgfi = p p U 2 p + Pg,p + PCR, P , (10) 

where P 9iP is the gas pressure at the point x = x p and Pcr, p 
is the pressure of accelerated particles at the same point 
(we use the symbol CR to mean cosmic rays, in the sense 
of accelerated particles). The mass and momentum escaping 
fluxes in the form of accelerated particles have reasonably 
been neglected. Note that at this point the equation for en- 
ergy conservation has not been used. 

Our basic assumption, similar to that used in (Eichler 
1984), is that the diffusion is p-dependent and more specifi- 
cally that the diffusion coefficient D(p) is an increasing func- 
tion of p. Therefore the typical distance that a particle with 
momentum p travels away from the shock is approximately 
Ax ~ D(p)/u v , larger for high energy particles than for 
lower energy particles*. As a consequence, at each given 
point x p only particles with momentum larger than p are 
able to affect appreciably the fluid. Strictly speaking the 
validity of the assumption depends on how strongly the dif- 
fusion coefficient depends on the momentum p (see section 
4). 

Since only particles with momentum > p can reach the 
point x — x p , we can write 



Pcr.p 



4-7T 



rpn 

J rj 



dpp v(p)f (p), 



(11) 



where v(p) is the velocity of particles with momentum p, 
Pmax is the maximum momentum achievable in the specific 
situation under investigation. 

From Eq. 10 we can see that there is a maximum dis- 
tance, corresponding to the propagation of particles with 
momentum p max such that at larger distances the fluid is 
unaffected by the accelerated particles and u p — u . 

The equation for momentum conservation is: 



dU 
dp 



1 _ 1 u-^+D 



Ml 



pou% dp 



(12) 



Using the definition of Pcr and multiplying by p, this equa- 
tion becomes 



dU_ 

' dp 



1 _ _Lr/-(7 S +i) 



' r p 4 v(p)Mp), 



3po«o 



(13) 



* For the cases of interest, D(p) increases with p faster than u p 
does, therefore Ax is a monotonically increasing function of p. 



where fo is known once U(p) is known. Eq. 13 is therefore 
an integral-differential nonlinear equation for U(p). The so- 
lution of this equation also provides the spectrum of the 
accelerated particles. 

The last missing piece is the connection between R S ub 
and Rtot, the two compression factors appearing in Eq. 8. 
The compression factor at the gas shock around x — can 
be written in terms of the Mach number Mi of the gas im- 
mediately upstream through the well known expression 



Rs 



ib = 



(7s + 



1)M? 



( 79 -l)M 1 2 + 2- 



(14) 



On the other hand, if the upstream gas evolution is adia- 
batic, then the Mach number at x = 0~ can be written in 
terms of the Mach number of the fluid at upstream infinity 
M as 



Ml ' = M ( ?(^) 73+1 =M 2 ( 
so that from the expression for R sn 

(7« + 1)^6 - (7» 



R 



sub 



Rtot J 
b we obtain 



R 



M 79+1 



1g + l 



(15) 



Now that an expression between R su b and Rtot has been 
found, Eq. 13 basically is an equation for R su b, with the 
boundary condition that U(p ma x) = 1. Finding the value 
of Rsub (and the corresponding value for R to t) such that 
U(pmax) = 1 also provides the whole function U(p) and, 
through Eq. 8, the distribution function fo(p)- If the reaction 
of the accelerated particles is small, the test particle solution 
is recovered. 



3 THE APPEARANCE OF MULTIPLE 
SOLUTIONS 

In the problem described in the previous section there are 
several independent parameters. While the Mach number of 
the shock and the maximum momentum of the particles are 
fixed by the physical conditions in the environment, the in- 
jection momentum and the acceleration efficiency are free 
parameters. The procedure to be followed to determine the 
solution was defined in (Blasi 2002): the basic problem is 
to find the value of Rsub (and therefore of Rtot) for which 
U(pmax) = 1. In Fig. 2 we plot U(p max ) as a function of R to t, 
for To = 10 s K, Pmax = 10 5 mc and pi„j = 10~ 2 mc in the left 
panel and p in j = 10~ 3 mc in the right panel (m here is the 
mass of protons). The parameter p was taken 10~ 4 in the 
left panel and 10~ 3 in the right panel. The different curves 
refer to different choices of the Mach number at upstream in- 
finity. The physical solutions are those corresponding to the 
intersection points with the horizontal line U(p ma x) = 1, 
so that multiple solutions occur for those values of the pa- 
rameters for which there is more than one intersection with 
U(p m ax) = 1. These solutions are all physically acceptable, 
as far as the conservation of mass, momentum and energy 
are concerned. 

It can be seen from both panels in Fig. 2 that for low 
values of the Mach number, only one solution is found. This 
solution may be significantly far from the quasi-linear so- 
lution. Indeed, for Mo = 10 the solution corresponds to 
Rtot ~ 8, instead of the usual R to t ~ 4 solution expected 
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Figure 2. Le/4 panel: U (p m ax) as a function of the total compression factor for T = 10 5 K, 77 = 10~ 4 , Pmax — 10 mc and £»i n - — 10 mc 
for the Mach numbers indicated. Right panel: same as in left panel but for r; = 10~ 3 and pi n j = 10~ 3 mc. 



in the linear regime. Lower values of the Mach number are 
required to fully recover the linear solution. 

When the Mach number is increased, there is a thresh- 
old value for which three solutions appear, one of which is 
the quasi-linear solution. For very large values of the Mach 
number the solution becomes one again, and it coincides 
with the quasi-linear solution. 

In Fig. 3 we show the appearance of the multiple so- 
lutions for the case To = 10 5 K, p ma x = 10 5 mc and p- n ~ = 

Trie with Mach number Mo — 100 {j>inj — 10 nxc and 
Mo — 100) in the left (right) panel. The curves here are ob- 
tained by changing the value of 77. The same comments we 
made for Fig. 2 apply here as well: low values of rj correspond 
to weakly modified shocks, while for increasingly larger ef- 
ficiencies multiple solutions appear. The solution becomes 
one again in the limit of large efficiencies, and it always cor- 
responds to strongly modified shocks. 

The problem of multiple solutions is not peculiar of the 
kinetic approaches to the non-linear theories of particle ac- 
celeration at shock waves. The same phenomenon was in 
fact found initially in two-fluid models (Drury & Volk 1980, 
1981), where however no information on the spectrum of 
the accelerated particles and on the injection efficiency was 
available. 



4 COMPARISON WITH AN ALTERNATIVE 
APPROACH 

Multiple solutions were also found by Malkov (1997) and 
Malkov et al. (2000), in the context of a semi-analytical 
kinetic approach. Aside from the technical differencies be- 
tween that method and the one proposed by Blasi (2002, 
2004), the main difference is in the fact that the former re- 
quires the knowledge of the exact expression for the diffusion 
coefficient as a function of the momentum of the particles, 
while the latter only requires that such diffusion coefficient is 
an increasing function of the particle momentum. While the 



first approach may provide us with an exact ' solution to the 
problem, the second is in fact more practical, in the sense of 
providing an approximate solution even in those cases, the 
majority, in which no detailed information on the diffusion 
properties of the fluid is available. The solution provided in 
(Blasi 2002, 2004) is particularly accurate when the diffu- 
sion coefficient is Bohm-like, D(p) = (1/3)t\lc, expected in 
the case of saturated self-generation of waves in the vicinity 
of the shock surface (Lagage & Cesarsky 1983). 

We will now discuss the quantitative comparison be- 
tween the results of Malkov (1997) and those of Blasi (2002, 
2004), by considering a single situation in which multiple 
solutions are predicted (in both approaches), and determin- 
ing the spectra and compression factors in both methods. 
We start with briefly summarizing the approach of Malkov 
(1997). The following flow potential is introduced there: 



dx'u(x'), 



(16) 



which is used as a new independent spatial variable to re- 
place x. Using the flow potential, it is possible to define an 
integral transformation of the flow profile as follows: 



U(p) = — 



110 



exp 



0+ 



q(p) 

~3D(p) 



du 

— da;, 
ax 



(17) 



where q(p) — —d\nfo/d\np is the spectral index of the par- 
ticle distribution function and D(p) is the diffusion coeffi- 
cient, which is assumed to be independent of the position. 
An integral equation for U (p) can be derived by applying 
Eq. 17 to the ^-derivative of the Euler equation (Malkov 
1997; Malkov et al. 2000): 

Rsub 



U(p) = 



Rt. 



- + 



Pin] J Pin. 



dp 



y/p 2 + (mc) 2 



1 + 



q(p)D(p) 
q(p)D(p) 



t Note however that even the approach of Malkov (1997) is based 
on several approximations: the solution is expanded to the first 
order and the contributions from gas pressure are ignored. 
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U (p i( 



U(p) 



■ exp 



RR S 



J Pin 



rilnp' 



where v is an injection parameter defined as 



47T 



C TPinjMPinj), 



(18) 



(19) 



3 po«o 

and related to the compression factor by the following equa- 
tion (Malkov 1997): 




exp 



Here R — Rtot/Rsub is the compression factor in the shock 
precursor. Eqs. 15, 18 and 20 form a closed system that can 
be solved numerically. 

Before showing the results it is worth noticing that the 
injection parameters r\ and v adopted in the two approaches 
are defined in two non equivalent ways. However, the relation 
between rj and v can easily be found by using Eqs. 8 and 19 
and can be written as: 



77 



Pi 



Rt. 



Rs 



1 



(21) 



One may notice that this relation contains the compression 
factors Rtot and R su b, which are what we are searching for. 
This fact implies that three solutions characterized by the 
same value of rj may correspond to three distinct values of 
v. 

In order to compare the results of the two different ap- 
proaches we consider a shock having Mach number Mo = 
150 and temperature at upstream infinity To = 10 4 K. We 
set the value of the injection and maximum momenta equal 
to 10~ 3 mc and 10 5 mc respectively and we adopt an effi- 
ciency 77 = 10 -4 . Using the approach proposed by Blasi 
(2002), we find three solutions, characterized by the values 
of the compression factor R ~ 15.3, 3.94, 1.05. The last so- 
lution is the quasi-linear one, in which the precursor is very 



weak. We adopt now these three values for the precursor 
compression factor to solve the system of equations given 
by Eq. 15, 18, and 20) for different choices of the diffusion 
coefficient. 

In Fig. 4 we plot the velocity profiles for the three solu- 
tions, as derived with the method of Blasi (2002, 2004) and 
detailed above (solid line). In the figure U(p) = u p /uo with 
Up defined through Eq. 5 for the method of Blasi (2002, 
2004). It is easy to show that U(p) is related to the U(p) 
through the relation U(p) + 1/Rtot = U(p). The dotted and 
dashed lines are the results obtained with the calculation of 
Malkov (1997) with a Bohm and Kolmogorov diffusion re- 
spectively. For Bohm diffusion the two approaches give very 
similar results. For Kolmogorov diffusion the differencies are 
larger, as expected. 

In Fig. 5 we plot the spectra of the accelerated particles, 
as obtained in this paper (solid line) and for a Bohm (dotted 
line) and Kolmogorov (dashed line) diffusion coefficient, as 
derived by carrying out the calculation of Malkov (1997). 

We recall that from the theoretical point of view the 
Bohm diffusion is in fact what should be expected in the 
proximity of a shock if the turbulence necessary for the ac- 
celeration is strong and generated by the same cosmic rays 
that are being accelerated (Lagage & Cesarsky 1983). In this 
perspective, we look at the results illustrated in this section 
as very encouraging in using the approach presented in Blasi 
(2002, 2004), since it is simple and at the same time accu- 
rate in reproducing the major physical aspects of particle 
acceleration at cosmic ray modified shocks. 



5 A RECIPE FOR INJECTION FROM THE 
THERMAL POOL 

The presence of multiple solutions is typical of many non- 
linear problems and should not be surprising from the math- 
ematical point of view. In terms of physical understanding 
however, multiple solutions may be disturbing. The typical 
situation that takes place in nature when multiple solutions 
appear in the description of other non linear systems is that 
(at least) one of the solutions is unstable and the system 
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Figure 4. Velocity profile upstream of the shock as derived in this 
paper (solid line) and with the approach of Malkov (1997) with 
a Bohm diffusion (dotted line) and for a Kolmogorov diffusion 
(dashed line). 
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Figure 5. Spectrum of the accelerated particles as derived in this 
paper (solid line) and with the approach of Malkov (1997) with 
a Bohm diffusion (dotted line) and for a Kolmogorov diffusion 
(dashed line). 



falls in a stable solution when perturbed. The stable solu- 
tions are the only ones that are physically meaningful. Some 
attempts to investigate the stability of cosmic ray modified 
shock waves have been made by Mond & Drury (1998) and 
Toptygin (1999), but all of them refer to the two-fluid mod- 



els. A step forward is being carried out by Blasi & Vietri (in 
preparation) in the context of kinetic models. 

In addition to the stability, another issue that enters 
the physical description of our problem is the identification 
of possible processes that determine some type of backre- 
action on the system. It may be expected that when some 
types of processes of self-regulation are included, the phe- 
nomenon of multiple solutions is reduced. In this section wc 
investigate the type of reaction that takes place when a self- 
consistent, though simple, recipe for the injection of particles 
from the thermal pool is adopted. This recipe is similar to 
that proposed by Kang, Jones & Gieseler (2002) in terms of 
the underlying physical interpretation of the injection, but 
probably simpler in its implementation. 

For non-relativistic shocks, the distribution of particles 
downstream is quasi-isotropic, so that the flux of particles 
crossing the shock surface from downstream to upstream can 
be written as 

-«d/ v (p) 



/oo /> — 1 

dp \ 



where v(p) is the velocity of particles with momentum p and 
Ud is the shock speed in the frame comoving with the down- 
stream fluid. The term Ud + v(p)n is the component along 
the direction perpendicular to the shock surface of the ve- 
locity of particles with momentum p moving in the direction 
Ijl. It follows that the flux of particles moving tangent to the 
shock surface (namely with p = —Ud/v{p)) is zero. We recall 
that, having in mind collisionless shocks, the typical thick- 
ness of the shock, A, is the collision length associated with 
the magnetic interactions that give rise to the formation of 
the discontinuity. Useless to say that these interactions are 
all but well known, and at present the best we can do is 
to attempt a phenomenological approach to take them into 
account, without having to deal with their detailed physi- 
cal understanding. It is however worth recalling that many 
attempts have been made to tackle the problem of injec- 
tion at a more fundamental level (Malkov 1998; Malkov & 
Volk 1995; Malkov & Volk 1998). Here, we consider the rea- 
sonable situation in which A oc r\ , where r^ h oc p t h is the 
Larmor radius of the particles in the downstream fluid that 
carry most of the thermal energy, namely those with mo- 
mentum 1.5 pth (pth = (2mfesT2) 1,/2 here is the momentum 
of the particles in the thermal peak of the maxwellian distri- 
bution in the downstream plasma, having temperature T2). 
We stress here the important point that the temperature 
of the downstream gas (and therefore p t h) is determined 
by the shock strength, which in the presence of accelerated 
particles, is affected by the pressure of the non-thermal com- 
ponent. In particular, the higher the efficiency of the shock 
as a particle accelerator, the weaker its efficiency in terms 
of heating of the background plasma (see section 8). 

For collisionless shocks, it is not clear whether the down- 
stream plasma can actually be thermalized and the distri- 
bution function be a maxwellian. On the other hand, it is 
generally assumed that this is the case, so that in the fol- 
lowing we consider the case in which the bulk of the back- 
ground plasma is thermal and has a maxwellian spectrum at 
temperature T given by the generalized Rankine-Hugoniot 
relations in the presence of accelerated particles (see section 
2). For modified shocks, the points discussed above apply to 
the so-called subshock, where the injection of particles from 
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the thermal pool is expected to take place. We recall that for 
strongly modified shocks the subshock is weak, and rather 
inefficient in the heating of the background plasma. 
From Eq. 22 we get: 

*=Jjf W , / tfc (p)^^, (23) 

where we assumed that the temperature downstream im- 
plies non-relativistic motion of the quasi-thermal particles 
(jp ~ mv(p)). In Eq. 23 we write the minimum momentum 
in terms of a parameter a, such that A = arf . With this 
formalism, the particles that can cross the shock surface are 
those that satisfy the condition: 

P> Pmin = 1.5 a pth- (24) 

The parameter a defines the thickness of the shock in units 
of the gyration radius of the bulk of the thermal particles. 
Our recipe for injection is pictorially illustrated in Fig. 6: 
thermal particles have a pathlength smaller than the shock 
thickness and cannot cross the shock surface, being advected 
downstream before the crossing occurs. Only particles with 
momentum sufficiently larger than the thermal momentum 
of the downstream particles can actually return upstream 
and be accelerated. 

In the following we will neglect the fluid speed Ud com- 
pared with v(p), which is a good approximation if the in- 
jected particles are sufficiently more energetic than the ther- 
mal particles. This is done only to make the interpretation 
of the result simpler, but there is no technical difficulty in 
keeping the dependence of the results on u d . 

We introduce an effective injection momentum pt n j = 
(,Pth defined by the equation: 

/>oo 

$= / dp4itp 2 f th (p)v(p), (25) 

which in terms of dimensionless quantities, with fth(p) = 
e -(v/p th ) 2 reads . 

poo poo 

/ dxx 3 e- x2 = 4 / dxx 3 e- x2 . (26) 

Jl.Sa J£ 

It is easy to show that £ « 2 for a = 1 (half a Larmor 
rotation of the particles with momentum 1.5p t h inside the 
thickness of the shock) and £ « 3.25 for a = 2 (one full 
Larmor rotation of the particles with momentum 1.5pth in- 
side the thickness of the shock. The fraction of particles at 
momentum £ times larger than the thermal one is ~ 5% for 
£, — 2 and ~ 10 -4 for £ = 3.25. The actual values of £ are ex- 
pected to be somewhat higher if the effect of advection with 
the downstream fluid is not neglected. The sharp decrease in 
the fraction of leaking particles that may take part in the ac- 
celeration process is due to the exponential behaviour of the 
maxwellian at large momenta. Although the fraction of par- 
ticles in the maxwellian that get accelerated only depends 
on the parameter £ which in turn is expected to keep the 
information about the microscopic structure of the shock, 
the absolute number of and energy carried by these par- 
ticles depend on the temperature of the downstream gas, 
which is an output of our calculation. This simple argument 
serves as an explanation of the physical reason why there 
is a nonlinear reaction on the system due to injection. If 
the parameter £ is assumed to be determined by the micro- 
physics of the shock, and if we adopt our simple recipe to 



describe such microphysics, then the shock thickness is eas- 
ily estimated once the temperature of the downstream gas is 
known, and the latter can be calculated from the modified 
Rankine-Hugoniot relations. The parameter r\ in Eq. 8 is no 
longer a free parameter, being related in a unique way to 
the parameter £ and to the physical conditions at the sub- 
shock. The condition that fixes ij is that the total number 
of particles in the non-thermal spectrum equals the number 
of particles in the maxwellian at momenta larger than pi n j . 
Due to the very strong dependence of the spectrum on the 
momentum for both the maxwellian and the power law at 
low momenta, the condition described above is very close to 
require the continuity of the distribution function, namely 
that fth(pinj) = /o(Pinj). In the following we adopt this 
condition for the calculations. This can be shown to imply 
the following expression for r\: 

r >=^(R™>>- 1 )t 3 e- e - (27) 

We recall that the compression factor at the subshock, R BU b, 
approaches unity when the shock becomes cosmic ray dom- 
inated. This makes evident how the backreaction discussed 
above works: when the shock becomes increasingly more 
modified, the efficiency r] tends to decrease, limiting the 
amount of energy that can be channelled in the non-thermal 
component. Although the recipe provided here is certainly 
far from representing the complexity of the reality of injec- 
tion of particles from the thermal pool, it may be considered 
as a useful attempt to include the main physical aspects of 
this phenomenon. 

6 SELF-CONSISTENT INJECTION AND 
MULTIPLE SOLUTIONS 

In this section we describe the role played by the injection 
recipe discussed above for the appearance of multiple solu- 
tions. It can be expected that the phenomenon is somewhat 
reduced because, as discussed in the previous section, the 
injection provides an efficient backreaction mechanism on 
the shock as a particle accelerator. Indeed we find that the 
appearance of multiple solutions is drastically reduced, and 
that the phenomenon still exists only in regions of the pa- 
rameter space which are very narrow and of limited physical 
interest. In the quantitative calculations we use the value 
£ = 3.5 for the injection parameter, as suggested by the 
simple estimate in Section 5 and as suggested also in the 
numerical work of Kang & Jones (1995). The dependence 
of the effect on the value of £ is discussed below. In Fig. 

7 we illustrate the dramatic change in the physical picture 
by plotting U(p m ax) as a function of R to t for £ = 3.5 and 
adopting the same values for the parameters as those used 
in obtaining Fig. 2. The efficiency r\ is now calculated ac- 
cording with the recipe described in the previous section. It 
can be seen very clearly that when the Mach number of the 
shock is changed, there is a single solution (compare with 
Fig. 2 where multiple solutions where found for the same 
values of the parameters, but without thermal leakage) . 

The appearance of multiple solutions can be investi- 
gated in the whole parameter space, in order to define the 
regions where the phenomenon appears, when it does. In 
Fig. 8 we highlight the regions where there are multiple so- 
lutions (dark regions) in a plane £ — log(p max ), for different 
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Figure 6. Graphic illustration of the structure of a collisionless shock and of the basic idea underlying our injection recipe. 
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uo = 5 x 10 8 cm s" 1 and p max = 10 3 , 10 4 , 10 5 , 10 7 mc from 
left to right. The line is continuous when there are no mul- 
tiple solutions and dashed when multiple solutions appear. 
The dashed regions are, as stressed above, rather narrow. For 
instance, for p ma x = W 4 mc there are multiple solutions only 
for 3.67 < £ < 3.7. Any small perturbation of the system 
that changes the values of £ at the percent level implies that 
the system shifts to one of the single solutions if it is sitting 
in the intermediate solution before the perturbation. The 
sharp transition between the strongly modified solution and 
the quasi-linear solution when £ is increased suggests that 
the intermediate solution may be unstable, though a formal 
demonstration cannot be provided here. In order to make 
sure that this is the case, a careful analysis of the stability 
is required, and will be presented in a forthcoming publica- 
tion (Blasi & Vietri, in preparation). On the other hand, a 
previous study, carried out in the context of the two-fluid 
models, showed that when multiple solutions are present, 
the solution with intermediate efficiency is in fact unstable 
to corrugations of the shock surface (Mond and Drury 1998). 



lot 



Figure 7. U(p ma x as a function of the total compression factor 
for T = 10 5 K, Pmax — 10 mc and £ — 3 for the Mach numbers 
indicated. 



ESCAPING FLUX OF ACCELERATED 
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values of the Mach number of the shock. In most cases the 
dark regions are very narrow and cover a region of values 
of ^ which is rather high (small efficiency). In Fig. 9 we 



It is rather remarkable that the kinetic model of Blasi (2002, 
2004) does not require explicitely the use of the equation for 
energy flux conservation. However, once the solution of the 
kinetic problem has been found, the equation for conserva- 
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Figure 8. Parameter space for multiple solutions. The dark regions illustrate the regions of parameters for which multiple solutions are 
still present. 



tion of the energy flux provides very useful information, as 
we show below. The equation can be written in the following 
form: 

-P2U 2 H —Pg.2U2 H ~P C ,2U2 = 

2 7s - 1 7c ~ 1 



1 3^ 

~Po«o + 

2 7s 



7s 



- -P 9 , o«o 



(28) 



where Fb is the flux of particles escaping at the maximum 
momentum from the upstream section of the fluid (Berezhko 
& Ellison, 1999). Notice that this term is usually neglected 
in the linear approach to particle acceleration at shock waves 
because the spectra are steep enough that, in most cases, we 
can neglect the flux of particles leaving the system at the 
maximum momentum. The fact that particles leave the sys- 
tem make the upstream fluid behave as a radiative fluid, and 
makes it more compressible. This is a crucial consequence of 
particle acceleration at modified shocks, and is shown here 
to be a natural consequence of energy conservation. 



In Eq. 28 we can divide all terms by (l/2)poMo an d 
calculate the normalized escaping flux: 



F'f. 



7 9 P< 



9.2 



R 2 tot Mj(7 9 -1) Rtot 7 9 - 1 Pou 2 



7c 



(29) 



Rtot 7c - 1 Pov% ' 
From momentum conservation at the subshock we also have: 

po«o Rtot Rtot jgM§ V Rtot ) 

so that the escaping flux only depends upon the environment 
parameters (for instance the Mach number at upstream in- 
finity) and the compression parameter R su b which is part of 
the solution. Note also that the adiabatic index for cosmic 
rays, 7 C , is here calculated self-consistently as: 



E c r ma!C dp4TYp 2 f (p)e(p) 



(31) 
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Figure 10. Escaping flux (F' E ), advected flux (F' adv ) and the sum of the two (F^ ot ) normalized to the incoming flux (l/2)po"o, as 
functions of the Mach number at upstream infinity Mq. Left panel: urj = 5 X 10 8 cm s -1 , p m ax = 10 6 mc and £ = 3.5. Right Panel: 
urj = 5 X 10 8 cm s _1 , Pmax = 10 2 mc and £ = 3.5. 
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Figure 9. Dependence of Rtot as a function of £ for Mo = 200, 
u = 5 X 10 8 cm s _1 and p ma x = 10 3 , 10 4 , 10 s , 10 7 mc from left 
to right. The sharpness of the transition suggests that the small 
perturbations of the parameters make the solution fall on one of 
the two sides. 



where E c is the energy density in the form of accelerated par- 
ticles and e(p) is the kinetic energy of a particle with momen- 
tum p. It can be easily seen that 7 C — > 4/3 when the energy 
budget is dominated by the particles with p ~ p ma x (namely 
for strongly modified shocks) and 7 C ^5/3 for weakly mod- 



ified shocks. In Eq. 29 the term F' adv = 



Rtot 7c -i 



Pa,2 



clearly the fraction of flux which is advected downstream 
with the fluid. 

In Fig. 10 we plot the escaping flux (F' E ), the advected 
flux {F' adv ) and the sum of the two (Fl ot ) normalized to the 
incoming flux (1/2)^0^0, as functions of the Mach number at 
upstream infinity Mo- Here we used uo = 5x 10 s cm s _1 , and 
£ = 3.5, while the maximum momentum has been chosen as 
Pmax = 10 6 mc in the left panel and p ma x = 10 2 mc in the 
right panel. Several comments are in order: 

1) At low Mach numbers the escaping flux is inessential, 
as one would expect for a weakly modified shock. We recall 
that the escaping flux is due to the particles with momentum 
Pmax leaving the system from upstream infinity. For a weakly 
modified shock at low Mach number the spectrum is steeper 
than E~ 2 , so that the energy carried by the highest energy 
particles is a small fraction of the total. 

2) At large Mach numbers the shock becomes increasingly 
more cosmic ray dominated, and for the cases at hand the 
total efficiency gets very close to unity, meaning that the 
shock behaves as an extremely efficient accelerator. At Mach 
numbers around 4 on the other hand the total efficiency is 
around 20% for p ma x = 10 6 mc and ~ 10% for p ma x = 
10 2 mc, dropping fast below Mach number 4. Clearly the 
efficiency would be higher in this region for lower values of 
the parameter £. 

3) Despite the fact that the total efficiency of the shock 
as a particle accelerator is close to unity at large Mach 
numbers, the fraction of the incoming energy which is actu- 
ally advected toward downstream infinity is only ~ 20% 
at Mo w 100 for p max = 10 6 mc. Most of the enegy 
flux in this case is in fact in the form of energy escaping 
from upstream infinity at the highest momentum p ma x- For 
Pmax = 10 mc the normalized advected flux roughly satu- 
rates at ~ 40% and is comparable with the escape flux at 
the same Mach number. For a distant observer these escap- 
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Figure 11. Temperature jump X2/T0 for p m ax /mc = 10 3 (th in 
solid line), Pmax/mc — 10° (dashed line), Pmax /mc = 10 7 (dot- 
ted line) and p m ax 

/mc = 5 X 10 1() (dash-dotted line). The thick 
solid line shows the jump for ordinary shocks. 
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Figure 12. Particle Spectra (thermal plus non-thermal) for Mo = 
10 (solid line), M = 50 (dashed line) and M = 100 (dotted line). 
The vertical dashed line represents the position of the thermal 
peak for ordinary shocks (it is almost independent of the Mach 
number for large Mach numbers). 



ing particles would have a spectrum close to a delta function 
around p max ■ 



8 SHOCK HEATING IN THE PRESENCE OF 
EFFICIENT PARTICLE ACCELERATION 

Energy conservation has the natural consequence that a 
smaller fraction of the kinetic energy of the fluid is con- 
verted into thermal energy of the downstream plasma in 
cosmic ray modified shocks, compared with the case of or- 
dinary shocks. The reduction of the heating at nonlinear 
shock waves is fully confirmed by our calculation in the con- 
text of the injection recipe introduced in section 5. In Fig. 11 
we plot the temperature jump between downstream infinity 
(at temperature T2) and upstream infinity (at temperature 
To). The thick solid line is the jump as predicted by the 
standard Rankine-Hugoniot relations without cosmic rays. 
The other lines represent the temperature jump at cosmic 
ray modified shocks with Pmax/mc — 10 3 (thin solid line), 
Pmax/mc = 10 5 (dashed line), Pmax/mc — 10 7 (dotted line) 
and Pmax/mc — 5 x 10 10 (dash-dotted line). 

Such a drastic reduction of the downstream tempera- 
ture is expected to reflect directly in the thermal emission 
of the downstream gas in those environments in which col- 
lisions are relevant. Note that for strongly modified shocks 
the compression factor between upstream infinity and down- 
stream are much larger than for ordinary shocks, so that the 
downstream turns out to be denser but colder than in the 
linear case. The missing energy ends up in the form of ac- 
celerated particles. 

The effect of suppression of the heating in cosmic ray 
modified shocks also appears in the spectra of the particles 
(thermal plus non-thermal) in the shock vicinity. In Fig. 



12 we show these spectra (including the maxwellian thermal 
bump) for uo = 5 x 10 s cm s _1 , £ = 3.5 and p max /mc = 10 5 . 
The vertical dashed line shows the position of the thermal 
peak as expected in the absence of accelerated particles. 
In fact this position should depend on the Mach number, 
but the dependence is very weak for large Mach numbers. 
The positions of the thermal peaks clearly show the effect of 
cooler downstream gases for modified shocks. At the same 
time, the effect is accompanied by increasingly more modi- 
fied spectra of accelerated particles, with most of the energy 
pushed toward the highest momenta. 



9 CONCLUSIONS 

The efficiency of the first order Fermi acceleration at shock 
fronts depends in a crucial way upon details of the mecha- 
nism that determines the injection of a small fraction of the 
particles from the thermal pool to the acceleration box. In re- 
ality the processes of formation of a collisionless shock wave, 
of plasma heating and particle acceleration are expected to 
be all parts of the same problem, though on different spatial 
scales. We hide our lack of knowledge of the microphysics of 
the shock structure in a simple recipe for injection, in which 
the particles that take part in the acceleration process are 
those that have momentum larger by a factor £ than the 
momentum of the thermal particles in the downstream fluid. 
This is motivated by the fact that for collisionless shocks the 
thickness of the shock surface is determined by the gyromo- 
tion of the bulk of the thermal particles. We estimated that 
£ ~ 2 — 4. This recipe implies that the fraction of particles 
that get accelerated is rather small (0.01 — 10~ 5 ). 

We implemented this recipe in a calculation of the non- 
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linear reaction of cosmic rays on the shock structure pro- 
posed by Blasi (2002, 2004). Similarly to other models, also 
this approach shows the appearance of multiple solutions, for 
a wide choice of the parameters of the problem. When the 
simple model of particle injection is used, this phenomenon 
is drastically reduced: the multiple solutions disappear for 
most of the parameter space, and when they appear they 
look as narrow strips in the parameter space, at the bound- 
ary between unmodified and modified shocks. We argued 
that this result suggests that the narrow regions may sig- 
nal the transition between two stable solutions, although 
this needs further confirmation through detailed analyses of 
the stability of the solutions. This interpretation seems to 
be supported in part by the calculations of Mond & Drury 
(1998), that showed that when three solutions are present, 
the intermediate one is indeed unstable for small corruga- 
tions of the shock front. This calculations was however per- 
formed in the context of a two-fluid model, while an inves- 
tigation of the stability for kinetic models is still lacking. 

We find that the phenomenology of the particle accel- 
eration at modified shocks is characterized by three main 
features: 

1) The modification of the shock increases with the Mach 
number of the fluid. For low Mach numbers the quasi-linear 
solution is recovered, but departures from it are evident al- 
ready at relatively low Mach numbers. The modification 
of the spectra manifests itself with a hardening at high 
momenta and a softening at low momenta. The p 4 /o(p) 
shows a characteristic dip at intermediate momenta, typi- 
cally around p/mc ~ 1 — 100 (for very large values of p m ax, 
the dip can be found at even larger momenta, which is of in- 
terest for the acceleration of ultra- high energy cosmic rays). 

2) The total efficiency for particle acceleration saturates at 
large Mach numbers at a number of order unity. However, 
as shown in Fig. 10, the largest fraction of the energy is 
not advected downstream but rather escapes from upstream 
infinity at the maximum momentum. This effect was also 
discussed in the context of the simple model by Berezhko & 
Ellison (1999). 

3) The high efficiency for particle acceleration reflects in a 
reduced ability of cosmic ray modified shocks in the heating 
of the background plasma. This effect is at the very basis 
of the backreaction introduced by the injection recipe on 
the acceleration process, and determines the saturation of 
the total efficiency for large Mach numbers. The heating 
suppression is shown in Fig. 11 and in Fig. 12. 
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